{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "fa89aa93-2169-4a16-9921-b91290a77052",
   "metadata": {},
   "source": [
    "## Prepare the surface fault rupture\n",
    "\n",
    "#### Reference\n",
    "+ Reitman, Nadine G, Richard W. Briggs, William D. Barnhart, Jessica A. Thompson Jobe, Christopher B. DuRoss, Alexandra E. Hatem, Ryan D. Gold, John D. Mejstrik, and Sinan Akçiz (2023) Preliminary fault rupture mapping of the 2023 M7.8 and M7.5 Türkiye Earthquakes. DOI: https://doi.org/10.5066/P985I7U2. Access date: 16 May 2023."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "ba949593-0a4c-4db0-bb68-c5eb6838b279",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Go to directory: /Users/yunjunz/data/archives/2023TurkeyEQ/USGS/simple_fault_2023-03-15\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import datetime\n",
    "import os\n",
    "import geopandas as gpd\n",
    "from matplotlib import pyplot as plt\n",
    "plt.rcParams.update({'font.size': 12})\n",
    "\n",
    "version = '2023-03-15'   # 2023-02-17, 2023-03-15\n",
    "work_dir = os.path.expanduser(f'~/data/archives/2023TurkeyEQ/USGS/simple_fault_{version}')\n",
    "os.chdir(work_dir)\n",
    "print('Go to directory:', work_dir)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7bb61862-1298-4d42-9f80-2def04722060",
   "metadata": {},
   "source": [
    "### Read shape file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "f2fbe001-96c6-4dc9-ad10-7a145dc65c59",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# read shape file\n",
    "suffix = datetime.date.fromisoformat(version).strftime('%Y-%-m-%d')\n",
    "shp_file = os.path.join(work_dir, f'simple_fault_{suffix}.shp')\n",
    "gdf = gpd.read_file(shp_file)\n",
    "line_strs = [geom for geom in gdf.geometry]\n",
    "line_ids = gdf.FaultTrace     # 1 is “confident” (solid lines) and 2 is “queried” (dashed lines)\n",
    "# remove lines with None\n",
    "line_ids = [x for x, y in zip(line_ids, line_strs) if y]\n",
    "line_strs = [y for y in line_strs if y]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2e4a97ce-8ea7-4b47-b87d-8f7fcc9db124",
   "metadata": {},
   "source": [
    "### Plot faults"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "32e2e9a9-7222-42a8-8a04-a4f7c4293716",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "save figure to file /Users/yunjunz/data/archives/2023TurkeyEQ/USGS/fault_2023-03-15.png\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhMAAAHRCAYAAADQazsGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB0hklEQVR4nO3dd3wUdf7H8dcmm957ICRA6L13pClSLRRFPRQsh4KeoqDY7gQ9f3bPs3A2REQURQFBpChFpClIkd5CIBAgvffs/P5YiMSEkMqmvJ+Pxz6S7MzOfHay2X1nvmVMhmEYiIiIiJSTna0LEBERkZpNYUJEREQqRGFCREREKkRhQkRERCpEYUJEREQqRGFCREREKkRhQkRERCpEYUJEREQqRGFCREREKkRhQqSG2LBhAyaTiZkzZ9q6lFqruh7jTz/9FJPJxKeffmrrUqpcXXqutYnChFSa+Ph4Pv74Y0aNGkXTpk1xcXHBy8uLvn37MmfOHCwWy2Ufu2XLFoYPH46vry+urq60b9+et956i/z8/CLrbt68mSeeeIJu3boRGBiIk5MTjRs35r777uPYsWPFbv+TTz7h5ptvpmnTpnh6euLm5karVq34+9//zuHDh0v1/CZOnIjJZCr1bcCAAaXabm1TmuM0derUYh+7Y8cO7r77bsLDw3FxccHT05MOHTowY8YMzp07V656Fi1axNChQwkMDMTBwQE/Pz9at27N+PHjmTdvXgWeae1ypd+ZrT/cGzVqRKNGjWxag1ye2dYFSO2xaNEiJk+eTHBwMIMGDSIsLIzz58+zePFi7rvvPn744Qe++eYbTCZTocd99913jBkzBmdnZ8aNG4evry/Lly/n0UcfZfPmzSxatKjQ+mPGjCE2NpZevXpxxx13YDab2bp1K3PmzGHhwoWsWbOG3r17F3rM559/ztmzZ+nRowfBwcHY2dmxf/9+5s6dy2effcbSpUsZNmxYic/v5ptvLvJmtmHDBn7++Wf69+9fJDzU9Te+m266iY4dOxa7rGfPnoV+NgyDJ598kldffRWz2czgwYO55ZZbyMnJYcuWLbz66qvMnj2bL7/8kpEjR5a6hkmTJvHRRx/h4uLCiBEjaNy4Menp6Rw/fpwlS5awYcMGJkyYULB+9+7dOXjwIP7+/uV6zrXBc889V+z9l/tdigBgiFSStWvXGkuXLjXy8vIK3X/27FkjNDTUAIxFixYVWpacnGz4+/sbjo6Oxvbt2wvuz8zMNHr16mUAxpdfflnoMS+99JIRFRVVZP8vvviiARht2rQpsiwzM7PYmtesWWMARsuWLUv9PC/13HPPGYDx3HPPlevxZbF+/fqrtq+KmDBhggEYc+fOLfVjZs6caQBGo0aNjH379hVZ/s033xjOzs6Gg4ODsXXr1lJt85dffjEAo0GDBsW+XtLS0ozvv/++1DXa0ty5c8t8TMsKMKrDR8LlnmvDhg2Nhg0b2qQmuTI1c0ilGTRoEDfddBP29vaF7g8ODuaBBx4ArP/JX2rRokXExcVx++2307Vr14L7nZ2d+fe//w3A7NmzCz3mySefpEGDBkX2P2PGDFxcXNi/fz9xcXGFljk7Oxdb8+DBg/H29ub48eOle5KlNHPmTEwmU5HnCxAZGYnJZGLixImF7r/YPBAREcFbb71Fu3btcHFxuWJzSVZWFmPHjsVkMvHggw8Wak768ssvGThwID4+Pjg7O9OqVSv+/e9/k52dXbBOYmIirq6uNGnSBOMyFxEeOXIkJpOJ33//vdTHoLROnDjBv//9bxwcHFi2bBlt2rQpss6YMWP4z3/+Q25ubsFr6Uo2b95c8NjiXi9ubm6MGDGi0H2X6zMxYMAATCYTubm5PP/88zRp0gRnZ2datGjBRx99VLDee++9R9u2bXFxcaFBgwbMnDmzSPPepb//Q4cOcfPNN+Pr64ubmxt9+/ZlzZo1pXp+F50+fZqHHnqI8PBwnJyc8PPz48Ybb2T79u1l2k5pREdH8/zzz9OnTx+Cg4NxdHSkfv363H777ezfv7/I+lfqg1KapouL2zh58iQnT54s1PTy178hsR01c8hV4ejoCICDg0Oh+9evXw/A0KFDizymX79+uLq6snXrVrKzs3FycipxHyaTCbPZ+pK++PVKNm3aRFJSEl26dCnV+lfDww8/zKZNmxgxYgTDhw8vEs4ulZiYyI033sjmzZt56aWXePLJJwuW3XvvvXzyySeEhoYyZswYvLy82LZtG//85z9Zu3Yta9aswcHBAR8fH2677Tbmzp3LTz/9xODBgwvtIyoqipUrV9KlS5cqOU5z584lLy+PW265hXbt2l12vfvuu4/nn3+ePXv2sHXrVnr16lXidgMCAgA4cuRIpdV622238euvvzJ8+HAcHBz45ptvmDRpEo6OjuzYsYMvvviCkSNHct1117F8+XJmzZqFi4sLM2bMKLKtEydO0KtXL9q2bcv999/P2bNn+eqrrxg2bBhffPEF48aNu2I9O3fu5PrrrychIYEhQ4YwevRo4uLiWLp0KX379mXJkiUMHz680p7/xo0befnllxk4cCBjxozBzc2No0eP8s0337Bs2TI2b95c6c0hjRo14rnnnuOtt94CKNTfRk0v1YfChFS5vLy8go5ufw0NFzs/NmvWrMjjzGYzjRs3Zv/+/URERNCqVasS97No0SJSU1Pp2bMn3t7exa7zzTffsG/fPjIzMzly5Ag//PAD/v7+vPfee+V4ZlVj165d7Nq1i8aNG5e43smTJxk2bBjHjh3js88+Y/z48QXLPv30Uz755BPGjh3L/PnzC52ZmTlzJrNmzeK9994reGN+8MEHmTt3Lh988EGRMPHRRx9hsVi4//77y/Q8li5dSmRkZLHLbrvtNlq2bAlYAx1QZL9/ZTabGTBgAF9++SW//PLLFcPE0KFD8fb2ZuXKldxwww3ceuutdOvWjebNm2NnV76TsqdOnWLfvn0Fr69p06bRsmVLHnnkEXx8fPjjjz8ICQkBYNasWTRt2pTXX3+dadOmFQm4GzduZPr06bz22msF9z300EP06tWLBx54gGHDhuHp6XnZWvLy8rj11ltJT09n48aN9O3bt2BZdHQ03bp145577iEyMvKyZ+aKU9xZhEaNGjFx4kQGDRrE+fPn8fDwKLR8586dXHPNNTz55JOsWrWq1PsqjUaNGjFz5syCDqDVbaSNXGDrdhap/aZNm2YAxrBhw4osa9asmQEYR48eLfaxvXv3NgBjy5YtJe4jIiLCCAgIMOzt7Y3Nmzdfdr1x48YVtA0DRrNmzYwdO3aU7Qld4nJ9Ji7ev379+iKPOXHihAEYEyZMKHT/xb4G//nPf4rd16V9Jnbt2mXUq1fP8PT0NH766aci63bs2NFwcHAwEhMTiyzLy8sz/Pz8jK5duxa6v1u3boaDg4Nx7ty5QuvWr1/f8PDwMNLS0oqt668uPo+SbkuWLClYv1WrVgZgrFy58orbfuKJJwzA+Mc//lGqWjZs2GA0bdq00L49PDyMYcOGGV9++aWRn59faP3L9Uvp37+/ARhr164tso+BAwcagDFnzpwiy+6++24DMCIjIwvuu/j79/LyMlJSUoo85uLx+/TTTwvuK64fwdKlSw3AePzxx4t97m+99ZYBlLpfSEm/r/79+1/x8SNHjjScnJyMnJycgvuu1M+nuH4Q6jNRM+nMhFSpt956izfeeIMWLVrw2WeflfnxxoU2/L+OALnU+fPnGTZsGLGxsbz99ttFRnJcauHChSxcuJCUlBT27dvHrFmz6N27Nx988EG1aX/t0aNHics3bdrEm2++iYeHBxs3bqRDhw6FlmdkZLBnzx78/f0LTg3/lZOTE4cOHSp035QpU7j77rv55JNPeOqppwBYvnw50dHRTJ48GTc3tzI9j7lz55bqmJbmd3zRxXWysrJKVUP//v05fPgwmzdv5ueff2bXrl1s3ryZlStXsnLlSj799FOWLVtW0Ax3JcU189SvX/+Ky06fPk3Dhg0LLevcuXOR//DB2j9j3rx57Nq1q9BIk7/aunUrYO2DUdx/60ePHgXg0KFDRfqGlMS4TL+Zi1asWMH777/Pjh07iIuLIy8vr9DyuLg46tWrV+r9Se2gMCFV5r///S+PPvoorVq1Yt26dcUOt/Py8gIgOTm52G2kpKQUWu+vzp07x6BBgzh8+DBvvfUW//jHP0pVm6enJ7179+b777+na9euTJ48meuuu67YjnpXW3BwcInLd+3aRWpqKr179y5oKrhUYmIihmEQGxvLrFmzSr3fcePGMW3aND7++GOefPJJTCYTH3zwAUCpOz2WR7169Th06BBRUVFXXPf06dPAn/0hSsPOzo5rrrmGa665BrB+WP74449MmDCB1atX87///Y9HHnmkVNsq7nV4sfmipGW5ublFlgUFBRW7j4u//8v9TVwUHx8PUGTo9F+lpaWVuLws3n777YImncGDBxMWFoarqysmk4mlS5eyZ8+eQp17pe7QaA6pEq+//jpTp06lbdu2bNiw4bIfkC1atACK7ySXl5fHiRMnMJvNhIeHF1keHR3NgAEDOHToEO+9916pPxAu5eDgwKBBg8jKymLbtm1lfvzlXGyT/+t/bQBJSUklPvZK/6E/9NBDPPDAA6xevZqbbrqJzMzMQssvfqh16tQJwzBKvF3KxcWFiRMnEhERwY8//khkZCRr1qyhZ8+etG/f/kpPudwutvX/9NNPJa6Xn59f0GG3Ih1BTSYT119/fcFoobVr15Z7WxVx/vz5Yu+/ODnX5QL0RReXf/fddyX+ji83b0RZ5eXl8dxzzxEcHMz+/fv56quveO2115g1axYzZ84sNhyV9HcAVw5MUnMoTEile+mll3j88cfp2LEj69evJzAw8LLrDho0CKDYTlsbN24kIyOD3r17FxnJERUVRf/+/Tly5Ajvv/8+U6ZMKXe9Z86cAUo/AqQ0fHx8AIr9b3vHjh0V2rbJZCr4b3r16tWMHDmS9PT0guXu7u60adOG/fv3k5CQUKZtT548ueCMRHk7XpbV3Xffjb29PYsXL+bAgQOXXe+TTz4hOjoaX1/fYkf/lNXFJoYrndavKjt37iQ1NbXI/ReHE3fq1KnEx1+c+OuXX36p9NqKExcXR1JSEr179y7SjJGWlsbOnTuLPKakv4Njx45dMVhfyt7evtgZcaV6UJiQSvXCCy/w9NNP06VLF9auXXvFmQTHjh2Lv78/CxcuLPQhm5WVxbPPPgtYP+AudfLkSfr3709ERARz5sxh0qRJJe4jPj6evXv3Frvs+++/Z8mSJbi7u9O/f//SPMVSudjv4eKwx4uioqJ4/vnnK2Ufb731FjNmzGDdunUMHTq00AfTY489Rk5ODvfcc0+xb9iJiYnFvvk3bdqUwYMHs2zZMj788EO8vb1LNUSxIho3bszTTz9Nbm4uN9xwQ7GBYunSpQVnnl555RVcXV2vuN1Vq1axePHiYpsY0tLSCvqT9OvXr2JPoJySk5OLvBZ27NjBggUL8PLyYtSoUSU+/qabbqJJkya89957/PDDD8Wus3XrVjIyMiql3sDAQFxdXdmxY0ehppPc3FweeeSRInO7ALRs2RJPT0++++47YmJiCu7PzMzk4YcfLtP+/fz8iI2NLXV/Gbm61GdCKs28efP417/+hb29Pddccw1vv/12kXUuDjG7yNPTk48++oixY8cyYMAAbrvtNnx9fVm2bBmHDx9m7NixRT7MBgwYQGRkJF26dOHkyZPFdj6bOHFiwWQ4UVFRdOrUic6dO9OmTRtCQkJISkpi9+7dbNu2DQcHBz7++OOC/6IqQ/fu3RkwYAAbNmyge/fuBUPqli9fzpAhQ0rVP6A0Xn75ZZycnHj++ee5/vrrWbVqFV5eXtxzzz38/vvvzJ49myZNmjBkyBDCwsJISEjgxIkTbNy4kbvvvpv333+/yDYnT57MmjVriIuL4+GHH8bFxaVctZU0NPSvr4OZM2eSnp7Om2++SYcOHRgyZAht2rQhNzeXLVu28OuvvwLwxBNPcN9995Vq/4cOHeLRRx/Fx8eHa665hmbNmmE2mzl9+jQrVqwgKSmJHj168NBDD5Xr+VVUv379+Pjjj/n111/p06dPwTwTFouFDz74oMRhoWBtolu8eDFDhgxhxIgR9O7dm44dO+Lq6kpUVBTbt28nIiKCs2fPlip8XYmdnR0PP/wwL7/8Mu3ateOmm24iJyeH9evXk5CQwMCBAwuaoS6t8bHHHmPmzJl06tSJUaNGkZeXx48//kj9+vULOqiWxrXXXsv27dsZNmwY11xzDY6OjnTo0IEbbrihws9NKsFVHTsitdrF4ZAl3S43xGzTpk3GsGHDDG9vb8PZ2dlo27at8eabbxaZmtswSh7CdvF26ZDMhIQE4+mnnzb69u1rBAcHGw4ODoarq6vRsmVL4/777zcOHDhQ4edc3NC3pKQkY9KkSUZAQIDh6OhotGnTxvjggw+uODT0xIkTxe6rpGF2//d//2cARufOnY34+PiC+5cvX26MGDHCCAgIMBwcHIygoCCjW7duxjPPPGMcPHiw2P3k5eUZ/v7+BmDs37+/tIeiyPMoz+tg27Ztxl133WU0bNjQcHJyKli/Xr16xo8//limOmJjY405c+YYt912m9GqVSvD29vbMJvNhr+/vzFgwADjvffeM7Kzsws95kpDQ0t6vsX93oobInzp7//gwYPGjTfeaHh7exsuLi5G7969jVWrVhXZTknTaZ8/f96YMWOG0aZNG8PFxcVwc3MzmjZtaowZM8aYP3++kZube8VjZRilm047NzfXeOONN4xWrVoZzs7ORlBQkDF+/HgjMjLyssfBYrEYr7zyihEeHm44ODgYoaGhxuOPP26kp6eXaWhoWlqa8cADDxghISGGvb19sX9DYjsmw7BRg6GIVEvHjx+nWbNm9O3bl40bN9q0ltTUVPr27cuBAwdYtGgRN998s03rqQyRkZE0btyYCRMm2PxKnCKVRX0mRKSQ1157DcMwbHb6/1IeHh58//33BAQEMG7cuEqfXVFEKof6TIgIJ0+eZP78+Rw9epT58+fTqVMnxo4da+uyAAgNDWXlypUsWbKEP/74g0GDBpV6kikRuToUJkSEEydO8M9//hM3NzeGDBnC//73v3Jfv6IqdOjQochMnyJSfajPhIiIiFRI9fnXQ0RERGqkWt/MYbFYiI6OxsPDo1QXEhIRERErwzBITU2lfv36JTZ91vowER0dTWhoqK3LEBERqbGioqJKvBBirQ8TF+ffj4qKuuKMciIiIvKnlJQUQkNDCz5LL6fWh4mLTRuenp4KEyIiIuVwpW4C6oApIiIiFaIwISIiIhWiMCEiIiIVojAhIiIiFaIwISIiIhWiMCEiIiIVojAhIiIiFaIwISIiIhWiMCEiIiIVojAhIiIiFaIwISIiIhWiMCEiIiIVojAhIiIiFaIwISIiIhWiMCEiIlXv+Dr4+i5IPWfrSqQKKEyIiEjVysmA7x+FA9/B1ndtXY1UAYUJERGpWhtfhcRI8AyB/jNsXY1UAYUJERGpOuf3w5Z3rN8Pfw2cPGxbj1QJhQkREakaabHw7d/BkgctR0LLEbauSKqI2dYFiIhI7ZOfcBL7z0dBwnFwC7CelZBaS2FCREQq1ZG92/FZPI4AIx68wuCupeBZ39ZlSRVSM4eIiFSan4/E8v2iOQQY8ZxxaAj3rga/JrYuS6qYzkyIiEil+Pb308z49g/yLCMJDPLgponTwTPI1mXJVaAwISIiFZKencerqw4xb+tJAG7uGMKtY1/B0ayT33WFwoSIiJTb5mNxzPj2D04nZgJwf/9wZgxpiZ2dycaVydWkMCEiImV2Kj6DN348zHe7owEI8XbhlTHt6dvM38aViS0oTIiISKnFpWXz7rpjLPj1JLn5BgB39mzIjGEtcXfSR0pdpd+8iIhcUURsGi98f4CtEfFk5VoA6Nc8gCeGtKBtiJeNqxNbU5gQEZFinYhLZ/X+c6zad47dUUkF97dv4MWTQ1vSu6maNMRKYUJERArsj05m1b5zrN5/jiPn0wotczLbMbpzA/5vVFtMJnWwlD8pTIiI1HF5+RZW7jvHJ5tPsOtUUsH9ZjsTvZr4cX2bYK5vHUSAu5NGaUixFCZEROqo5IxcFm4/xbwtkUQnZwHgaG/HoJaBDGkbxKAWQXi5Oti4SqkJFCZEROqY5Ixc/vPTEb7eEUVGTj4Afm6OjO/ZkPE9GxLg4WTjCqWmUZgQEakjDMNg2Z5oXvj+AHFpOQC0DPbgnr6NubFDfZwd7G1codRUChMiInVAVEIGzyzdx8YjsQA0DXTnXyNbc00zf3WmlApTmBARqcVy8y18/MsJ/rv2CFm5FhzNdjw0sCn39w/HyawzEVI5FCZERGqpXacSeWrxXg6dSwWgV7gfL45qS3iAu40rk9pGYUJEpJZJycrl9dWHmb/tJIYB3q4OPDO8FWO7NFCThlQJhQkRkVrCMAxW7TvHzOX7OZ+SDcDoziE8M7wVfu4aoSFVR2FC6qx8i0FqVi6pWXmkXPhqveWSlp1HVm4+OXkWcvIsZOdbyMs3cHGwx9XJHjdHM66O9rg6mgt+dnGwx8nBDmfzn19dnexxsLez9VOVWi4jJ4+dJ5P48JeIgg6WjfxceXFUO/poymu5ChQmpNayWAzScvLYdyaZfWeSOR6TTkxqFrFp2cSkZBOfnkO+xajSGkwmCPJwJtTXhVAfV0J9rbcmAW60DPbExVEd4KRiEtNzuP4/G4lNs56JMNuZeKB/Ex4a1FRDPWuD7FSIOwKxRyDhOGQmQU6a9f7cDPBvAc2HQMM+YHa0WZkmwzCq9t3UxlJSUvDy8iI5ORlPT09blyNVyGIxOBabxvbIBH6LSGDV/nNk51mu+DhnBzs8nB3wcDbj4eyAp7MZdyczzg72ONrb4Wi23sx2JrLzLKRn55GRk096Th4Z2flk5OaRnp1PVq71lp1nISs3nyvlFDsTNPJ3o3U9T1rX96R1PU+aBLhTz8sZs85mSCmcScrk+eX7Wb3/PAADWgQw68Y2NPRzs3FlUm6JJ+HYT3B8HUTvgpQzpXucowc0GQgthkGH263/yVSC0n6G6syE1Fg5eRb2nE5ie2QCOyIT+f1kIsmZuUXWq+flTMdQb1oGexLs5USAhxOBHs4EeDjh4+qIo7nyP7gNwyA33yA5M5fTiRlEJWYSlZDB6cQMTiVkcOR8GrGp2UTEphMRm873f5wteKzZzkQDHxdCfV0Ju+QW6utKeIAbro76s63rcvMtvLf+GP/bcJzsPAt2JnhyWEsm9Wti69LqhLx8C6cSMsjJt5BvMQj1dcXTuZzTjqfHw8nNELnJGiDijxZdxy0QAlqAX1Nw8wdHd3ByB3tHOPUrHF0D6TFwcBkkn4aOd1TsCZaD3pWkxjAMgxNx6Ww8EssvR+PYGhFfMBXwRS4O9oQHuLE/OgWAN2/twOjODa56rSaTCUeziQAPa3jpFOZTZJ2Y1CwOnk3lQHQKB86mcCA6maiETHLyLUTGZxAZn1HMdqGBjwvNAj1oFuhO00B3mgV50DTQHXcn/TnXBdFJmfzjy138fjIRgB6NffnXDa1pU9/LxpXVDTsiE5i+aE+hv0+TCVoEeTAoIIURLTxoE+pv/aC3Mxc9Q2DJh3N/QOSFABF7sPBykz2Edoem10LDvhDYElyKvn8U6HwXWCxwdhccWQ1eoZX4bEtPzRxSraVn57HpWBwbDsey8UgsZ5IyCy33c3Oke2NfujbypVsjH+LSsnlwwS4yc/O5tWsDXh3bwUaVl0++xeB8ShanEqxnMKIufD0Zn8HJ+HQSM4qeebkoxNvFGi4C3WkW5E6b+l60CPZQB9BaZP3hGB77ajeJGbl4OJt5cVQ7bmhfT8M9q1hUQgZbjsex8WgcP+w9i2FYL8dutjORfsk/NP9z+A/D7LeXfQcBraBRH2jcD8IHgHP1CYZq5pAa60xSJusOnuengzFsjYgn55J+D472dnRt5MM1zQLo19yfVsGeBZdEzs230O3Fn8jMzad/8wBm3djWVk+h3OztTNT3dqG+tws9w/2KLI9Py+ZYTBpHYtI4dj6VozFpHI2xNpmcScrkTFImP1/ozQ/gaLajTX1POjTwpn0DLzqEetPYz02Xka5h8vItvPnjEWZvOA5A2xBPZt/RhTA/VxtXVjvFpGax+VgcW47Fs+V4fJF/YgAshkHuX7pkJRnuxJj88HM2YW/kQn5e8TvwbQyN+lo7TTbsbW26qOEUJsTmLBaDfdHJ/HTgPD8ejOHg2ZRCy0N9Xbi2ZRD9mvvTM9zvsn0GTsZnkJSRi6ujPR/d1bVK+kLYmp+7E37uTvT4S9BIysjh2IVgcfR8GkfOp/LH6SRSsvLYdSqJXaeSCtb1cDLTur4nbUO8aBviSbsQLxr7u2OvgFEtJabn8I8vd7HpWBwAd/ZsyDMjWmmkRiWLSc3i+z1nWbnvLDtOJnLpOXuznYmOod70buJHSlYen26JJDffAAxCvF0Y3DqIQS0Dad9gMN6uthtRYUsKE2ITefkWfotMYNW+c6zef65ggh2wjnLoHObDta2CuLZVIM0C3Ut1GjciNg2w9imojUGiJN6ujnRtZG3uucgwDE7GZ7DndBJ7opL543QS+6KTSc3O49cTCfx6IqFgXVdHe1rXswaMNvU9aRbkQWN/N7xcytmpTCrMMAy2HI9nxrd/cDoxExcHe14Z254bO9S3dWm1RmpWLqv3n+e73WfYfCyu0AisdiFe9G7qR+8m/nRt6IPbhT5JyZm53NEjDLA2dYT5uqqZiTKGid27d/PMM8+wd+9eYmNjcXFxoUWLFjz44IOMHz++YD3DMPj44495//33OXr0KA4ODrRt25YnnniCESNGlGpfP/30E//85z/Zs2cPrq6ujBw5kldffZXAwMCyPUOpNiwW65vjir3RrNl/nvj0nIJlbo729GsewLWtghjYIqBcs/WFB7hjMsGR82n8fjKRLg1L6LRUB5hMJhr5u9HI342bOoYA1hB3PDadvRfm3th3JpkDZ1PIyMlnx8lEdlzo1HeRn5sjjf3drLcAN0J9XAnxcaGBtwv+7k5qLqkCGTl5LNl1hs+2nOTwees1NcJ8Xfnwri60DFa/r4rKybOw4XAM3+2O5qeD5wsNH+8U5s0N7eszrF0w9bxcin28l4uDQnYxyhQmkpKSCA0N5fbbbyckJIT09HQWLFjAnXfeSWRkJM8++ywAzz33HC+88AIPPPAAL7/8MllZWbzzzjuMHDmSb7/9ltGjR5e4n59//plhw4YxYsQIvvvuO2JiYpgxYwbXXnstO3bswMlJ08LWJBk5eSzeeYa5m09wPDa94H5vVwcGtwpiWLtg+jT1r/AVDJsGujOyfX2W77G+SdT1MFEcs70dLYI9aBHswdgu1lEu+RaDE3Fp7DuTwr4zyeyPTiEiLo3zFyb2ik/PKRIywNp/pZ63M/W9XAjxsfbzaODtQgMfFxr4uFLP21mdP8vgVHwGn22N5OsdUaRkWdvaXRzsGdMlhOnXt6izp88rS3ZePu+tP868LZGFhpCHB7hxc8cQbupYX/NzVECljObo2bMn0dHRnDp1CoAGDRrQuHFjfvnll4J1srKyCA4Opn///nz33Xclbq979+6kp6ezZ88ezGZr3tmyZQt9+vRh9uzZTJ48udS1aTTH1ZeVm8/pxAx2Homk8e8v8UziCI5keQPg7mTmxo71Gd62Hj3CfSv9w+b/fjjIhxsjuL9fOE8Nb1Wp265r0rLziIxLJyIunROx6ZyIS+N0YibRSZmcS8kq1aRcwZ7ONPBxpYGPNXCEXOhY2shfb9pgPVu36Vgc87ZEsu5wTEE7fUM/V+7s2ZBbuobqv+BKsPNUIk988wfHYqxNoYEeTtzYoT43dwqhTX1PNVOU4KqO5vD39ycmJqbgZwcHB7y8Cg9tcXZ2LriV5MyZM2zfvp2XXnqpIEgA9O7dm+bNm7NkyZIyhQmpOoZhkJadR0xqNkfPp7LhcCybjsVxOjGTTqajvOP4Dg1McczMP8FTvi8yoXdjbunaAI/yTu5SCqlZ1v849OZQce5O5gudNIsOU8vNt3A+JYszidYRJNEXRpKcTvzza06ehejkLKKTs/gtsvDjG/m50r2xL90a+dKjsR+hvi515ndmsRjsOJnID3vP8sPes8Sk/tlfqF/zACb2bsiA5oFqQqoE0UmZ/Peno3z9exSGAf7uTsy8sTXD2tZTh+NKVq4wYbFYsFgsJCYmsmjRIlavXs27775bsPyRRx5h+vTpzJkzh9GjR5OVlcVrr71GcnIyDz/8cInb3rdvHwDt27cvsqx9+/Zs3ry5xMdnZ2eTnf3nH2dKSkoJa0t53PrBVvafSSazmCmjTViYbP8908xfYzZZSHQKweGaWazrPbDK/3izcvNZue8cAD3Cfa+wtlSEg73dhTMOxQ9NtFgM4tKzOZ2YeeGWQXRSJsdi0tgRmVgwKdfXO04D1v8UO4R60zHUm85hPvRo7FurPkxLChAezmbGdG7Anb0a0iTA3YZV1h6J6TnM3nCMeVtPFgwtH905hH+NbK3moipSrjAxZcoUPvjgAwAcHR15++23uf/++wuWT506FRcXFx588EHuu+8+AHx9fVm+fDl9+vQpcdvx8fEF6/+Vr69vwfLLeemll5g1a1aZno+UTXZufqGJWtwc7Qn1daV3Yy/uT3iVoJPfWxe0GY3PDW/R7SpNwLJk1xmSMnIJ8XahX7OAq7JPKZ6dnYlAD2cCPZzp/JfZP1OyctkRmcBvJxLZHpnAH6eTiEnN5scD5/nxgPUaE438XJnYuxFju4bW2Jk9DcNgf3QK3/x+utgAcX3rYEa0r5z+QmKVnp3HJ5tO8OHGCFKzrf1Oujf2ZcbQFnRpqH8wqlK5/kqffvpp7rvvPmJiYli+fDkPPfQQ6enpTJ8+HYC5c+fyyCOP8NBDDzFs2DBycnL47LPPuOmmm1i8eDFDhgy54j4ud8rzSqdCn3rqKR577LGCn1NSUggNtc30orXV27d3wjCswwndnc3WeR9yM+Hru+DkGusUsiPegM4TKu1iM1eSlp3Hmz8eAeDuPo10CrMa83R2YFDLIAa1DAKsZ5T2nUlmd1QSu6OS2Hgklsj4DGYuP8Aba45wS9dQ+jT1o0mAO6G+rtX+d2sYBqv2nePtdccKzZmiAFF1cvIsfPnbKd5Zd5S4NOsosdb1PHliaAv6Nw+oM01otlSuMBEWFkZYmHWc7fDhwwHrh/iECRMwm80FZyRef/31gscMGzaMAQMG8MADD3DixInLbtvPzzoZT3FnIBISEoo9Y3EpJycnjfaoYsX2eF79jPViM2ZnuHU+NL/+qtb0wc/HiU3NppGfK3f1anRV9y0V4+xgX2iOjPTsPBbvPM3cLZFExKbzyeYTfLLZ+p7haLYj3N+NQE9nLBaDPIsFiwXcnOxp4ONKq3qeNPZ3Y3tkAtsi4om4ZPQQWGcY7RjmzaAWgfRvEYB/OYYgl2TLsTheWXWIPaeTC+q9vnUQozqF0LeZAkRly87LZ9nuaP679iinE62zVDb0c2Xa9S0Y2a5erWoqq+4q5fxh9+7def/994mIiMAwDDIzM+nWrVuR9bp27crPP/9MWloa7u7Ftw22bWudAnnv3r0FQeWivXv3FiyXauTkFtgxx/r9bQug6XVXdfcWi8Enm6wfNk8Oa1nnJqyqbdyczNzZqxF/69GQjUdjWbLrDEfOpxERm0Z2noVD51I5dC613Ns/k5TJij/OYjJB+wbe9G8eQOt6njQPcqehn1uZznxcnBhs83HrNRs2H7P+E+TqaM9914Rzb5/GeLlqNEZli0nJ4vNfT/HFrycLzkQEeDjxyLXNGNctVEOSbaBSwsT69euxs7MjPDyczExrOty2bRsTJkwoWMcwDLZt24aPjw9ubpcfFhYSEkL37t35/PPPmT59Ovb29gXbO3z4MFOnTq2MkqWy5GbBsn9Yv+9811UPEgD5hlHQh6NXeM2f416s7OxMDGgRyIAW1onq8i0GZxIzORabSkJ6LmY7E3Z2JuxNJlKycolKyODnI7HEp+XQpZEPPcP9aBfihfmScJCWncfmY3GsPxzDvjMp7IlKYk9UUsFyR7MdTQIuXCwt0J1ATye8XR3xcXXEx9UBb2c7LCZ7tkXEs/lYHJuPFb5ug4O9ib/1aMiDA5sS4KEzpJVtd1QSn24+wYq9Zy9MZ20dfjyxTyMm9GqEi6PO/NhKmcLEpEmT8PT0pHv37gQFBREXF8eiRYv46quvePzxxwkIsHZ6Gz16NB9++CFOTk4MHz6c7Oxs5s2bx+bNm3nhhRcKtV+ZzWb69+/P2rVrC+575ZVXGDx4MLfccgtTpkwhJiaGJ598krZt23L33XdX0lOXSmHJs16wJicdBr9g62qkFrO3MxHm51rixa2eGNryitvpGe7HtOtbcD4li58Px7I1Ip6jMakci0kjK9fCwbMpRa4PUxIHO4NOQWb6NPVndI+WhPprREZlys7LZ9W+c3y6JbLQNWa6NvTh7j6Nub5NkM5EVANlmrRq7ty5zJ07l4MHD5KUlIS7uzsdOnTgvvvuKzSddlZWFu+++y7z58/nxIkTODg40Lx5cx566CHuuOOOQmHCZDLRv39/NmzYUGhfP/74I//617/YvXt3wXTar732Wpmn09akVVdJZhK4eNtk1/kWg+bPriTfYrDi4b60qV99Lt8rNcfFMx9HzqdyJCaViNh04tOySczIJSkjh8SMXFIycwCD1qaT9LHbTx+7fXSzO4yr6cJIDXsn8GkIPo3Bp5H16pBNBkFAC1s+tRolN9/C3jPJbD0ez7aIeHZEJpKZaz3z6Ghvxw0d6jOxdyPaNdDf+dVQ2s/QSpkBszpTmKgbpiz4nR/2nmNUpxD+M66jrcuRWio/M5Xc+BM4p56ExEhIOAGJJ6xfk6OsZ+r+6oa3ocuEovdLgbi0bJbtjmbj0Vi2n0goNPQcIMjTib/1aMjt3cPUfHSVXdUZMEVsbXL/pvyw9xzL9kTz2ODmhPpe/lS4SHnZu3hg36A9UHRSPfLzIOX0hYAR+WfICFan8eLk5Vv4+UgsX++IYu3BGPIumQHP29WBHo196RXuR68m/jQPKt2Vg8V2FCakVmjXwItrmvnzy9E45mw6wcwb29i6JKlr7M3Wpg2fRraupFrKtxgcjUll16kkdp1KZMPh2EITeXUI9eaG9vXo1cSPVsGeGtZZwyhMSK1xT9/G/HI0jqW7z/DU8JYa0y9iQwnpOfx+MpFdpxLZdSqJP04nFWm+8HVzZHSnEG7pGkqLYA8bVSqVQWFCao1+zQII9nTmXEoWPx2IYUT7erYuSaROOZecxap9Z1m57xzbIxOKXLvHzdG+4BosXRv50LdpgOaFqSUUJqTWsLczMaZLCO+tP863O08rTIhcBVEJGay8ECAuHboJ0DTQnc5h3nQK86FTmDfNAj2q/XToUj4KE1KrXNMsgPfWH+dEXPqVVxaRMsvJs7D3TBKbj8Wzev859kcXnpOjS0MfhrUNZkibYHWErkMUJqRW0X89IpUrN9/CH6eT2RZRdN4HADsT9Gjsx7B21gAR5Olsw2rFVhQmpFa5GCUstXv6FJEqlZiew/d/RLPmwHl+P5lIxl86Tvq4OtAz3I/+zQMY3DoIv0q+YJrUPAoTUqt4uzoCkHDh4j8iUjq5+RY2HI7l299Ps/bQ+YJrX4A1PPRo7EfPcF96NfGnWaC7hm5KIQoTUqsEe1lPsaZm55GenYebk17iIiXZH53MN7+fZtnuaOLT/wzhret5MqpTCNc096d5oIfCg5RI77RSq7g7mfFwNpOalceBsyl0a+Rr65JEqh2LxWD94Rg+3BjBrycSCu73d3diVKf6jO7cgFb1dPkBKT2FCal1rm8dzLc7T/PRxgiFCZFLZOXms3TXGT76JYLjsdYRT2Y7E0PaBjO2cwOuaeaPWVfglHJQmJBaZ/KAcL7deZqtEfEkZeQU9KMQqYuy8/LZdyaFjUdiWfDrSeIu9CfycDJzR88wJvZuRD0vFxtXKTWdwoTUOk0DPZj9t870aeqPl4uDrcsRuari07LZczqJ7ZGJ/B6ZyO7TSeTkWQqW1/dy5p6+jRnXLRQPZ/19SOVQmJBaaXg7zX4ptZdhGEQnZ3EsJq3gdjwmjWOxaSSkFx3J5OfmSJeGPoxoX4/h7erhoKYMqWQKEyIiNUBieg4/HTzP6v3n2Xo8rshFsy4VHuBG14Y+dG3kS9eGPjT2d9MlvKVKKUyIiFRT8WnZLNsTzer959gemUj+JVfOMtuZaOTvRtMAd5oGutMsyJ0mAdabi6OumCtXl8KEiEg1cyYpk482RrBw+ymycv/s79CqnidD2gRxXasgWgR7qLlCqg2FCRGRauJYTBrv/3ycpbvOkHfhLES7EC9u6lhfF86Sak1hQkTExvaeTmb2hmOs2n+Oi5eV6dPUjykDmtK7iZ/6O0i1pzAhImIDhmGwLSKB2RuO8cvRuIL7r28dxJSBTekY6m274kTKSGFCROQqOxCdwvPf72dbhHUqa3s7Ezd1qM8DA5rQPMjDxtWJlJ3ChIjIVXLoXArvrD3Gir1nAXA02zGuayiT+oVXTX8IwwA1kchVoDAhIlLF4g5v5b+/pjD/QF7BfSPb12PG0JalDhGGYZCSlUdieg6JGTnYHV8L8cdo6ZKIU1o0ZCVDdirkpFm/ZqdZv3d0B1dfcPMHV7/ib95h4NcEHN2q6hBILacwISJSRSwWg8/W7WToL3cwjWz+sHuSBm368tCgple8KmdEbBrbIhI4cj6VQ+dSOHK+8OyWqxxn0tIu6spF5KRab0knr7yuR31rqPBr+ufNv5k1bNhr6m25PIUJEZEqkJWbz2Nf7WL04ekE2ydwxr4BL98zllYN6xdaLzkzl83H4jidmEFCei7HY9M4eDaF04mZxW7XzdEeb1dHDtCFZFMjmjZrjV9IE3DxAScP65kIJw/rzcEVcjMgIx7S46xfM+IhIw4yEi7cHwuJkdbvU6Ott8hfCu/Uvzk8tL2KjpTUBgoTIiKVLDE9h79/toN2p7/gOodd5Ns5Uu/eL/ALCOLg2RT2RCWxZNcZjpxPJTEjt9ht2JmgZ7gfbUO8aB7kQctgj7/Mbjmo9AX5NbnyOhkJkBABcUch/tiF23HrV5/Gpd+X1EkKEyIilehkfDoT527HLX4vzzguAOBzr0l8+Fk80cmrCuaRuFSTADfahXjh5eJAI383WgZ70rqeJ16uV7FpwdXXemvQtfD9Fou174VICRQmREQqwmKBHXOg/Th2x1qYMXc1w7LXco/TasxYyDdMvHG2PSlYmy08nc20CPagf/MABrYMpIG369UNDWVlZwfOJffvEFGYEBEpr+xUWDoZDi4nedtnJMSbWcFuzA7W62lk4sw+74E817kJoWHhhAe44efmqBktpdZRmBARKY/447DwbxB7kHzMeCX8waALGSG/QU/su07ApfVNdHN0o5ttKxWpcgoTIiJldfQnjG/vwZSVTAw+PJD9CI+YF2Nfrw09x0zFHNTS1hWKXFUKEyIipWUYsPktjJ9mYcLgd0szJudMxa9eQ9xumkDXRr62rlDEJhQmRERKKfOX93BZNxMT8EXeQF4338cjN7blbz3CMNvb2bo8EZtRmBARKYXIo/sIXjcLgJdzbyO+4xTWDGuJv7uTjSsTsT2FCRGRK9hwOIbtX77P4+Sww64dg//+El3UpCFSQGFCROQyDMPgk82RvLjiABZjGDn1mjFlzFB8GihIiFxKYUJEpBjnkrN4eeVBlu6OBuDWrg14/OZhOJrVN0LkrxQmREQukZyRy/9+Ps7czSfIzrNgZ4Knh7fi3r6NNdmUyGUoTIiIYL3K56dbIpm9/hgpWXkAdGnow9PDW9KloZo1REqiMCEidd76wzE8vXgvZ5OzAGge5M4TQ1pybatAnY0QKQWFCRGps5Izc/n39wdY9PtpAEK8XXh0cHNGdQrB3k4hQqS0FCZEpE5afziGp77dy7mULEwmuLt3Yx4f0gIXR3tblyZS4yhMiEid8tezEY38XHntlg5007wRIuWmMCEidYbORohUDYUJEan1jsek8t6G4yzeeQbQ2QiRyqYwISK12kcbj/PSykNYDHQ2QqSKKEyISK1ksRj87+fjvL76MAbgYG/HJxO7ck2zAFuXJlLrKEyISK2TlJHDY1/vYd2hGAB6N/HjvTs64+PmaOPKRGonhQkRqVX+OJ3ElAU7OZ2YiZPZjhduasut3UJtXZZIraYwISK1gmEYfPHbKWYtO0BOvoUwX1f+N74zbep72bo0kVpPYUJEaryMnDyeXbKPxbusozUGtw7i9Vs64OXiYOPKROoGhQkRqdGOx6Yx5fOdHD6fir2diSeGtGBSv3BdU0PkKlKYEJEa64e9Z3nimz9Iy84jwMOJd2/vRI9wP1uXJVLnKEyISI2Tm2/hpR8O8cnmEwD0aOzLO3d0ItDD2caVidRNChMiUmPkWwx+2HuW/649yrGYNAAe6N+E6dc3x2xvZ+PqROouhQkRqfYuhoi31x7l6IUQ4e3qwGtjOzC4dZCNqxMRhQkRqbaKCxGezmb+fk04E/o0wtNZozVEqgOFCRGpdhQiRGoWhQkRqTYUIkRqJoUJEbE5hQiRmk1hQkRsRiFCpHZQmBCRq+5yIeK+a8KZqBAhUuMoTIjIVWOxGKxQiBCpdRQmROSq+DUinhd/OMgfp5MBhQiR2kRhQkSq1PHYNF5eeYgfD5wHwM3Rnr/3C+eevo0VIkRqCYUJEakS8WnZ/HftURb8eop8i4G9nYnbu4fyyLXNCfBwsnV5IlKJFCZEpFIlZ+QyZ/MJPtl0grTsPACuaxXIk8Na0jTQw8bViUhVUJgQkUpxMUTM3XSC1Ashom2IJ08Pb0XvJv42rk5EqpLChIhUSHJGLnM2RTB3c2RBiGgZ7MEj1zZjSJtg7OxMNq5QRKqawoSIlItChIhcpDAhIqWWlZvPxiOxrNx3jh8PnC/oE6EQIVK3KUyISImOxaSyPTKRX47GsuFwLBk5+QXLFCJEBBQmROQyTsVnsOC3k3y88QT5hlFwf30vZ4a2rcewdsF0CfNRiBCRsoWJ3bt388wzz7B3715iY2NxcXGhRYsWPPjgg4wfP75gPZPp8m8uLVq04NChQyXuZ8CAAfz8889F7h8yZAirVq0qS8kiUkoWi8Gx2DR+PhzLir1n2R2VVGj5mC4h3NmzER0aeJX4Ny4idU+ZwkRSUhKhoaHcfvvthISEkJ6ezoIFC7jzzjuJjIzk2WefBWDr1q1FHvvrr78ydepURo0aVap9hYeHs2DBgkL3eXt7l6VcESlBenYee6KS+P1kIjtOJrLrVCIpWXkFy00muKapP5uOxWEx4PHrWxLs5WzDikWkujIZxiXnL8upZ8+eREdHc+rUqcuuc/fddzNv3jyOHDlC06ZNS9zegAEDiIuLY9++fRUtjZSUFLy8vEhOTsbT07PC2xOpiQzD4ExSJr+fTCy4HTybguUvf/0uDvZ0bujN0DbBDGkbjKO9HR2f/xGAI/8ehqPZzgbVi4itlPYztFL6TPj7+xMTE3PZ5ampqSxatIj+/ftfMUiISMXl5FnYH53M7ycT2XnKGh7Op2QXWS/E24XODX3o2tCHLg19aBnsgdn+z8Bw7MKVPT2czQoSInJZ5QoTFosFi8VCYmIiixYtYvXq1bz77ruXXX/hwoWkp6dz3333lXofx48fx9fXl5SUFBo2bMhtt93Gs88+i4uLS4mPy87OJjv7zzfNlJSUUu9TpCaLSsjgy99OsT0ygT9OJ5OdZym03Gxnok2IF13CrMGhc0Nv6nmV/PcUlZABgL+7rqUhIpdXrjAxZcoUPvjgAwAcHR15++23uf/++y+7/pw5c/D29mbMmDGl2n7fvn0ZN24cLVu2JDMzk5UrV/Lqq6+yadMm1q9fj53d5f9Deumll5g1a1bZnpBIDXYmKZN31x1j0Y4o8i5pt/BxdbgQGnzoEuZD+wbeuDjal2nbX++IAqBvU02HLSKXV64+E6dOnSImJoaYmBiWL1/Ohx9+yCuvvML06dOLrLt//37atm3Lgw8+WOLZiyt54403mD59OosXLy6xE2dxZyZCQ0PVZ0JqneikTN5bf4yvd0SRm2/9M+7b1J8bO9ana0MfGvu7VWjUxbnkLPq8so58i8Hqqf1oEayLdInUNVXaZyIsLIywsDAAhg8fDsBTTz3FhAkTCAgIKLTunDlzAMrUxFGc8ePHM336dLZt21ZimHBycsLJSadkpfY6m5zJ7PXH+Wp7FDn51qaM3k38eHRwc7o18q20/Xy4MYJ8i0H3xr4KEiJSokrpgNm9e3fef/99IiIiCoWJnJwc5s+fT5cuXejYsWNl7KrEJg6R2ux8Shaz1x/jy9/+DBE9w3159Lrm9Aj3q9R9HT6XyrytkQA8NFCdpkWkZJUSJi72YwgPDy90/7Jly4iLi+P555+v8D7mzZsHWIehitQlx2JSmbMpkm93nibnQqfK7o2tIaJXk8oNEWAdRvqv7/aRbzEY0iaIfs0DrvwgEanTyhQmJk2ahKenJ927dycoKIi4uDgWLVrEV199xeOPP15sE4eLiwt33HHH5Qswm+nfvz9r164F4JdffuHFF19k1KhRhIeHk5WVxcqVK/nwww8ZNGgQN9xwQzmepkjNYhgGG4/GMWfTCTYeiS24v1sjn4IQUVWzUC7bE82vJxJwMtvxz5Gtq2QfIlK7lClM9OrVi7lz5zJv3jySkpJwd3enQ4cOzJ8/v9B02gBRUVGsWbOG8ePH4+Xlddlt5ufnk5//54WD6tWrh729PS+88AJxcXGYTCaaNWvG888/z7Rp09TMIbVaVm4+S3ad4ZNNJzh6YY4Hkwmubx3EvX3D6dbIp0qnsjYMg7d+OgrAgwOb0sDHtcr2JSK1R6XMgFmdaQZMqQliUrOYv/UkC349RUJ6DgBujvbc2i2Uu3s3Jszv6nyo7z2dzA3vbsLZwY7fnx2Mm5OuBShSl13VGTBFpHz2RyczZ9MJlu+JLhjeGeLtwt19GnFrt1A8nR2uaj3f740G4NqWQQoSIlJqercQucosFoO1h2KYsymCbREJBfd3aejDvX0bc33roEJTWl8t+RaDFX+cBWBk+3pXff8iUnMpTIhcJenZeXzz+2nmbj5BZLx1mmp7OxPD29Xj3r6N6RjqbdP6Fu2I4nRiJp7OZga0CLRpLSJSsyhMiFSx5IxcZm84xpe/nSq4xLens5nbe4QxoVcj6nuXfH2MqyE5M5fXVh8G4JHrmpd52m0RqdsUJkSqUERsGvfO28GJuHQAGvu7cXefRozp3KBa9Ul4e+1R4tNzaBrozl29Gtq6HBGpYarPu5lILbP1eDwPfP47yZm5hHi7MOvGNgxqGYidXdUN7SyPrcfjmbclEoB/jmyNgw36a4hIzaYwIVIFvt4RxdOL95JnMegY6s1Hd3UlwKP6XTPmWEwq98/fQZ7F4KaO9emv2S5FpBwUJkQqkcVi8Orqw7z/83HAOiri9Vs64OxQPfsgPLdsPylZeXRp6MMrY9rbuhwRqaEUJkQqSVZuPo99vZsf9p4D4OFBTZl6XfNq16xxUXp2Hr+dsA5NfW1s+2obeESk+lOYEKkECek5/P2zHfx+MhFHezteHtOO0Z0b2LqsEv12IoHcfIMGPi409nezdTkiUoMpTIhU0Mn4dCbO3c6JuHQ8nc18eFdXelbyJcGrwi9H4wDo29S/Sq/3ISK1n8KESDkZhsGhP37l4NLX6J9TnxzvUXx6dzeaBXnYurRS2XzsQpho5m/jSkSkplOYECmnl1cdYtnGo3ztsIcWLtFMmfIagZ7Oti6rVOLSsjl8PhWA3k0UJkSkYhQmRMohIjaND3+OwMCX7yy9eaCpgbmGBAmAbRHxALQM9sDXzdHG1YhITafZaUTKyDAM7v50Owbga5fOffY/YPaoWdey2HrcGiZ6Nan+fTtEpPpTmBApo/lbT3LywoW6Hm96GmdTLtjXrP/uL4YJNXGISGVQmBApg4jYNF5YcQAAbxcHbmmQauOKyi46KZOIuHTsTNC9sa+tyxGRWkBhQqSUsnLzmbJgJ7n5BgCPXNcMM3kXltacoZUbj8QC0DHUGy8XBxtXIyK1gcKESCn9e8UBDp2znonwdDEzrlsoJJ6wLvQKsWFlZfPzhTDRv3nN6uchItWXwoRIKXz/RzSfbztV8PP4Hg1xdTTDeWuTB/4tbFRZ2cSnZReEiX7N1V9CRCqHwoTIFZyMT+epb/cCcPEyGzd1DIHk05BwHEx20KCrDSssvXfXHyMjJ592IV50DPW2dTkiUksoTIiU4HxKFhPnbic1O4/G/q5YDGga6E7zIHeI2GBdqX5ncPG2ZZmlciYpkwUXzq7MGNpSU2iLSKVRmBC5jNjUbO74aBsn4tJp4ONCiLcrAMPb1bN+EF8ME00G2q7IMgjycOLfN7dldOcQTaEtIpVKYUKkGPFp1iBxPDad+l7OzL+3O1GJ1rklRrSrB7lZcGytdeXwAbYrtAzM9nbc2i2UN2/taOtSRKSW0XTaIn+RnJnLnXN+42hMGsGeznw5qScN/dxYP20Ae88kW5s4dn4GmQng2QBCe9i6ZBERm9KZCZFLZObkc++n2zlwNgV/d0e++HsPGvq5AWBnZ6JDqDem/FzY8o71AT0ng73mahCRuk1hQuSCnDwL93/+OztOJuLpbGb+vT0ID3AvvJJhwA/TIf4oOHtB57tsU6yISDWiMCEC5FsMpn61i41HYnFxsGfu3d1pVc+z6Iq/vg875wEmGP0ROBezjohIHaMwIXWeYRg8vXgvP+w9h6O9HR/e1YUuDX2KrnhsLax+2vr99S9A8yFXt1ARkWpKYULqvJdXHuKrHVHYmeDt2ztyTbOA4ldc/yIYFuj4N+j10NUtUkSkGlOYkDrtfxuO88HGCABeHt2eoW3rFb+iYUDsEev3fR4BTfgkIlJAYULqrGV7onll1SEAnh7eklu7hV5+5fQ4yEkFTODd8OoUKCJSQyhMSJ10Ii6dp779A4BJ/cKZ1K9JyQ8ouDpoA3BwruLqRERqFk1aJXVObr6Fh77YSXpOPt0b+/LEkFJc8dOzPgz5P0DNGyIif6UwIXXOnE0n2B+dgo+rA2/f1gmzfSlO0Hk1gF4PVn1xIiI1kJo5pE45nZjBf386CsAzI1oT7KUmCxGRilKYkDrDMAz+uXQfmbn59Gjsy5jOIbYuSUSkVlCYkDpj8c4zrD8ci6PZjhdHtbVeRlxERCpMYULqhJiULGYt3w/A1Oua0TTQw8YViYjUHgoTUutZLAZPLt5LSlYe7UK8mHRNuK1LEhGpVRQmpNZ7ZfUh1h2KwdHejtduaV+60RsiIlJqeleVWu2r7af44GfrdNmvjm1Py2Bd5VNEpLIpTEitteV4HM8s2QfAw9c24+ZOGr0hIlIVFCakVjp6PpXJn+8kz2JwQ4f6PHpdM1uXJCJSaylMSK3z+8lEbvlgK8mZuXQK8+a1se01DFREpAppOm2pVXZHJfG3j7eRlWuhY6g3n0zohrODva3LEhGp1RQmpNYwDIPnlu0nK9fCNc38+eDOLrg66iUuIlLV1Mwhtcb3f5xlT1QSro72vHFrBwUJEZGrRGFCaoXsvHxeXX0IgPv7NSHQQxfwEhG5WhQmpFb4dHMkUQmZBHo48fd+jW1djohInaIwITVeREwar6yynpWYPKCJmjdERK4yhQmp0fZHJ3PDe5uxGGBvZ2JM5wa2LklEpM7Rv3BSI+XlW/h252lmLT9ARk4+AI9e1wxPFwcbVyYiUvcoTEiNYhgGK/ed4/U1h4mITS+4v4GPC5P6NbFhZSIidZfChNQYsanZ3PfZDvZEJRW6v2mgO++P74KjWa12IiK2oDAhNcas5fvZE5WEvZ2JfIsBwOhOIfx7VFt1uhQRsSG9A0uNsPlYHN//cRaAfIuB2c7Ecze2YXyPMF13Q0TExhQmpNpLSMtm0mc7Cn72dXNk9t860zPcz4ZViYjIRQoTUq3l5lsY9b8tpF8YsdEiyIOPJ3Ql1NfVxpWJiMhFChNSbRmGwT++2MXJ+AwAbu8Wyj9vaK3+ESIi1YzelaXaenvtUVbtPwdA2/qe/N/oduofISJSDWksnVRLPx04z39+OgqAg72J/43voiAhIlJNKUxItZOckcvTS/YW/Dz1uubqIyEiUo0pTEi18/z3B4hJzQbA392Re/roKqAiItWZwoRUK5uOxvHtztMFPz9ybTNcHO1tWJGIiFyJwoRUK19uP1XwfaivC+O6hdmwGhERKQ2FCak2dkclseLCLJcA9/drouttiIjUABoaKjZnGAbrD8fwnx+tozfMdibyLAadw3xsXJmIiJSGwoTY1LGYVJ78di87TiYC4GS2IzvPgr2diSaBbjauTkRESkPnkMVmsvPyuefTHew4mYizgx339GnM62M7ABDk4YSTWR0vRURqAp2ZEJv5ZFMkpxIyCPJ0YumDfajn5cL6QzEAeLk62rg6EREpLZ2ZEJuITc3mvfXHAHhiSEvqebkAcOR8KgDhAWriEBGpKRQmxCa8XR14bHBzrmnmz6hOIYC1I+aaA+cBaBXsYcvyRESkDNTMITbhYG/HPX0bc3efRgXX3Fiy6wy/n0zExcGe0Z0b2LhCEREprTKdmdi9ezcjRowgLCwMFxcXfH196dWrF59//nmh9Uwm02VvLVu2LNW+fvrpJ3r16oWrqyv+/v5MnDiRmJiYspQrNcDFIJGZk88rqw4B8I9rm1Lf28WWZYmISBmU6cxEUlISoaGh3H777YSEhJCens6CBQu48847iYyM5NlnnwVg69atRR7766+/MnXqVEaNGnXF/fz8888MGzaMESNG8N133xETE8OMGTO49tpr2bFjB05OTmUpW2qATzaf4HxKNiHeLroWh4hIDWMyDMOo6EZ69uxJdHQ0p06duuw6d999N/PmzePIkSM0bdq0xO11796d9PR09uzZg9lszTtbtmyhT58+zJ49m8mTJ5e6tpSUFLy8vEhOTsbT07PUj5OrJy4tmwGvbSAtO4+3xnXk5gt9KERExLZK+xlaKR0w/f39Cz70i5OamsqiRYvo37//FYPEmTNn2L59O3feeWehbfbu3ZvmzZuzZMmSyihZqpH/++Egadl5tA3x5MYO9W1djoiIlFG5OmBaLBYsFguJiYksWrSI1atX8+677152/YULF5Kens599913xW3v27cPgPbt2xdZ1r59ezZv3lzi47Ozs8nOzi74OSUl5Yr7FNvZFhHP4p1nMJnghZvaYmdnsnVJIiJSRuU6MzFlyhQcHBwIDAzk0Ucf5e233+b++++/7Ppz5szB29ubMWPGXHHb8fHxAPj6+hZZ5uvrW7D8cl566SW8vLwKbqGhoVfcp9hGdFImUxfuBuCO7mF00rU4RERqpHKFiaeffprt27ezYsUK7rnnHh566CFef/31Ytfdv38/v/76K3/7299wdnYu9T4u9vIv7f0XPfXUUyQnJxfcoqKiSr1PuXqSM3OZOPc3zqVk0SzQnRnDSjfKR0REqp9yNXOEhYURFhYGwPDhwwHrh/iECRMICAgotO6cOXMAStXEAeDn5wdQ7BmIhISEYs9YXMrJyUmjPWqAGd/8wZHzaQR6OPHpPd3xdHawdUkiIlJOldIBs3v37uTl5REREVHo/pycHObPn0+XLl3o2LFjqbbVtm1bAPbu3Vtk2d69ewuWS80VlZDBqv3nAPhkYjdCNKeEiEiNVilhYv369djZ2REeHl7o/mXLlhEXF8e9995b6m2FhITQvXt3Pv/8c/Lz8wvu37ZtG4cPH2b06NGVUbLY0IcbraGzb1N/2oZ42bgaERGpqDI1c0yaNAlPT0+6d+9OUFAQcXFxLFq0iK+++orHH3+82CYOFxcX7rjjjssXYDbTv39/1q5dW3DfK6+8wuDBg7nllluYMmUKMTExPPnkk7Rt25a77767jE9RqpMzSZks3G6dj+TBgSUPExYRkZqhTGGiV69ezJ07l3nz5pGUlIS7uzsdOnRg/vz5jB8/vtC6UVFRrFmzhvHjx+Pldfn/PvPz8wudgQAYMGAAP/zwA//617+44YYbcHV1ZeTIkbz22mvqD1HDvbvuKLn5Br2b+NGriZ+tyxERkUpQKTNgVmeaAbP6OBWfwaA3NpBnMfjmgV50bVRyZ1oREbGtqzoDpkhpvL3uKHkWg37NAxQkRERqEYUJuSpOxKWzeOdpAB4b3NzG1YiISGVSmJCr4p11R7EYMLBFAB1DvW1djoiIVCKFCalyEbFpLN11BoCp1+mshIhIbaMwIVXunXXHsBhwbctAOuishIhIraMwIVXqbHIm3+3WWQkRkdpMYUKq1O8nE7EY0C7Ei3YNNNuliEhtpDAhVWrjkVgAOoV527YQERGpMgoTUmWycvNZudd6Qa/h7erZuBoREakqChNSZbZGxJOanUc9L2e6a5IqEZFaS2FCqszmo3EA9G8egJ2dycbViIhIVVGYkCqz6Zg1TPRp6m/jSkREpCopTEiViEvL5tC5VAB66+qgIiK1msKEVIlfjlpHcbSq54mfuy4bLyJSmylMSJVYd8gaJga1DLBxJSIiUtUUJqTS5eRZ+PlwDACDWgbZuBoREalqChNS6d7/+TgpWXn4uzvpCqEiInWAwoRUqoNnU3hn3VEA/jmyFfYaEioiUuspTEilyc23MH3RHnLzDQa3DuLGDvVtXZKIiFwFChNSaTYeiWV/dApeLg68OKotJpPOSoiI1AUKE1Jpdp5KBODaloEEejjbuBoREblaFCakUmTl5vPV9igA+rfQcFARkbpEYUIqxdJdZ4hLy6G+l7OuECoiUscoTEiFWSwGH286AcA9fRvjYK+XlYhIXaJ3famw9YfPcywmDQ8nM+O6hdq6HBERucoUJqTCHv16DwA3dKiHh7ODjasREZGrTWFCKmTr8ThSMvMA+FvPhjauRkREbEFhQipke6R1OGh4gBtt6nvZuBoREbEFhQkpt8T0HD7aGAHAI9c2s3E1IiJiKwoTUm7vrj9GanYerep5ckN7TZ0tIlJXKUxIuUQlZDB/60kAZgxtgZ0u6CUiUmcpTEi5vLHmMDn5Fvo09aN/c814KSJSlylMSJkdOpfC0t3RADw1rJUu6CUiUscpTEiZ/efHIwCMaFePtiEawSEiUtcpTEiZ/BoRz+r957EzwdTrNIJDREQUJqQM8i0Gs5YfAOC27mE0C/KwcUUiIlIdKExIqX29I4oDZ1PwcDYzbXBzW5cjIiLVhMKElEpadh6vrz4MwNTrmuPn7mTjikREpLpQmJBSWbLrDPHpOTTyc+WuXroGh4iI/ElhQq7IYjGYvzUSgAm9G+Fgr5eNiIj8SZ8KckWfbD7BkfNpuDnaM7pzA1uXIyIi1YzChJToyPlUXr3QV+KZEa3xcnGwcUUiIlLdKEzIZeXkWXj0q93k5FkY2CKA27uH2rokERGphhQm5LLeXXeU/dEpeLs68MqY9po2W0REiqUwIcXaEZnAu+uPAfDize0I9HS2cUUiIlJdKUxIEYfOpvC3j3/FYsBNHeszon09W5ckIiLVmNnWBUj1YRgGC349xazl+8nNN7AzwTPDW9m6LBERqeYUJgSAM0mZPLV4LxuPxAIQ5uvK1OuaqXlDRESuSGGijjMMgy9+O8VLPxwiLTsPR7MdM4a25O7ejbCzU4dLERG5MoWJOiw2NZupX+1i87F4ALo09OHVse1pEuBu48pERKQmUZioo07FZ3DnJ79yMj4DZwc7Hh/Skom9G2GvsxEiIlJGChN10IHoFO765Dfi0rJp4OPCvHu662yEiIiUm8JEHZOSlcuEudYg0TLYg8/u6a5OliIiUiEKE3XMm2uOEJuaTbi/G1/d30vX2hARkQrTpFV1yL4zyXx24VLiz9/UVkFCREQqhcJEHWHJz+f5pTuxGDCifT36NvO3dUkiIlJLKEzUEbvXLeS9mAn83XEN/xzR2tbliIhILaI+E3VE/X0fEmBKpodfNsFe6nApIiKVR2cm6oJTvxKcvJtsw8y2gFttXY2IiNQyChO1nWHAL68DsDS/L9kugTYuSEREahuFidpuz0I4uoZ8kz0f5o/AwV6/chERqVz6ZKnNEiPhh8cB2FjvPo4bITia9SsXEZHKpU+W2sowYMlkyEmFsF785H8HAM4O+pWLiEjl0idLbXVoBZzaAg6uMOoDUrMNANydNIBHREQql8JEbWQYsOEl6/c9J4NPQ9Kz8wDwcFaYEBGRyqUwURud3g7n94HZBXr/A4DYtGwAPJw1hbaIiFQuhYnaaOc869c2N4OLDylZuew7kwxAh1Bvm5UlIiK1k8JEbZOdCvsWW7/vPAGAXyMSsBjQ2N+NEG8XGxYnIiK1kRrQaxtHdxi/GI6shLCeAHy0MQKALg29bViYiIjUVgoTtY3JBA17WW/At7+f5rfIBADMdjoRJSIilU+fLrWYYRjMWr4fgK6NfJgxtIWNKxIRkdpIZyZqsfj0HFKy8jCZ4Iv7emr2SxERqRL6dKnFTidkABDo7qQgISIiVUafMLXYhsOxACRn5mKxGDauRkREaqsyhYndu3czYsQIwsLCcHFxwdfXl169evH5558XWTc3N5c333yTdu3a4eLigre3N71792bLli1X3M+AAQMwmUxFbkOHDi1LuXVaVm4+3/8RDUCgpzN2diYbVyQiIrVVmfpMJCUlERoayu23305ISAjp6eksWLCAO++8k8jISJ599lkA8vPzGTVqFJs2beKJJ56gd+/epKen8/vvv5Oenl6qfYWHh7NgwYJC93l7e5el3DorPi2bSfN/51hsOmY7mP23TrYuSUREajGTYRgVPv/ds2dPoqOjOXXqFABvvfUW06ZNY/PmzfTs2bPM2xswYABxcXHs27evoqWRkpKCl5cXycnJeHp6Vnh71V1adh43vLOJE3HpeDqbeX98F3o39bd1WSIiUgOV9jO0UvpM+Pv7Yzb/eZLjv//9L/369StXkJCK+WhjBCfi0qnn5cziKX0UJEREpMqVK0xYLBby8vKIjY1l9uzZrF69mhkzZgAQFRVFZGQk7dq14+mnnyYoKAiz2UybNm2YN29eqfdx/PhxfH19MZvNNGnShGeeeYbMzMwrPi47O5uUlJRCt7oiJjWLj36xznb57IjWNA10t3FFIiJSF5RrnokpU6bwwQcfAODo6Mjbb7/N/fffD8CZM2cAmDdvHg0aNODdd9/Fy8uLjz76iIkTJ5KTk8Pf//73Erfft29fxo0bR8uWLcnMzGTlypW8+uqrbNq0ifXr12NXwkyOL730ErNmzSrP06rRDMPghe8PkpGTT4dQb4a3C7Z1SSIiUkeUq8/EqVOniImJISYmhuXLl/Phhx/yyiuvMH36dLZs2UKfPn1wdHTkyJEjNGzYELB+2HXt2pWYmBiioqLKXOgbb7zB9OnTWbx4MaNGjbrsetnZ2WRnZxf8nJKSQmhoaK3vM/HZ1kj+9d1+7O1MLHqgF53DfGxdkoiI1HBV2mciLCyMrl27Mnz4cP73v/8xadIknnrqKWJjY/Hz8wOgZcuWBUECwGQyMWTIEE6fPk1MTEyZ9zl+/HgAtm3bVuJ6Tk5OeHp6FrrVdrtOJfLC9wcAeGpYSwUJERG5qiqlA2b37t3Jy8sjIiKCJk2a4OrqWux6F0+ClNRMcSUVeWxtlJqVyz++3EVuvsGwtsHc27exrUsSEZE6plI+mS/2YwgPD8dsNnPTTTdx8OBBIiMjC9YxDINVq1bRpEkT/P3LPsLgYudNjRAp7PnlBzidmEkDHxdeHdsek0mTU4mIyNVVpg6YkyZNwtPTk+7duxMUFERcXByLFi3iq6++4vHHHycgIACAF154gZUrVzJ06FBmzpyJp6cnH3/8MXv27OHrr78uXIDZTP/+/Vm7di0Av/zyCy+++CKjRo0iPDycrKwsVq5cyYcffsigQYO44YYbKump13wpmbl8+/spwI43bumAh7ODrUsSEZE6qExholevXsydO5d58+aRlJSEu7s7HTp0YP78+QV9GgCaNGnCL7/8wpNPPsmkSZPIzc2lY8eOLFu2jJEjRxbaZn5+Pvn5+QU/16tXD3t7e1544QXi4uIwmUw0a9aM559/nmnTpqmZ44LcvHyWfPg8PzguYWzOTDqGedu6JBERqaMqZQbM6qxWzoCZFEXGN5NxPf0LAPM97uPOaW/YuCgREaltSvsZWq55JsSGTmyEr+7ENSuJTMORT5zvYtLDL9u6KhERqcPUZlCT/P4pzB8FWUnstoQzIvclbrz/BRwc1FdCRERsR2cmaoqNr8G6fwNwLGgo407eRpuwQEJ9ix+GKyIicrUoTNQEkZth3YsAZPWZwS1bupBNHqM7N7BxYSIiImrmqP4yk2DJ/YABHcfz7/QbSMzMo1mgO7d1C7V1dSIiIgoT1d4P0yE5Cnwas6vtDD7fdgqAWTe1wWyvX5+IiNiePo2qsz8Wwd5FYLKH0R/x6Y54AMZ0bkDvJmWfRVRERKQqKExUV0mnYMVj1u/7P0FKQEdW7TsHwITeDUt4oIiIyNWlMFEdZSbBwjsgOwUadMfSZxqPLtxNdp6F5kHutAvxsnWFIiIiBTSaozo5vQMiNsCR1XBuL7gFwJiPmLM1irWHrJdtv751sC7mJSIi1YrCRHWx/WP4YQYYedafnTzhzqVE5PnzyqqNBau1b6CzEiIiUr2omaM6OLcPVkz7M0gA/G0RBLfly99OkWexXj7F0d5E54Y+NipSRESkeDozUR1EbbN+DWoL5/eBVyiE9SQnz8K3O88UrHZ//yb4uzvZqEgRqalyc3MLXZ1ZxN7evlIvxaAwUR2c/cP61dXP+rVeBwAc7E2M7hTCx5tO4O/uyAP9m9ioQBGpiVJSUoiLiyM7O9vWpUg15OTkhL+/f6VcUVthojo4dyFMZKdYv9bvBEBadh6Ld1nPTDw6uDluTvp1iUjppKSkcObMGdzd3fH398fBwUGdtwUAwzDIzc0lOTmZM2esnzEVDRT6dLK1/Fw4f8D6fdxR69fwgQB89MsJEtJzCPd349aumjpbREovLi4Od3d3GjRooBAhRbi4uODh4cHp06eJi4urcJhQB0xbizsC+dlgdoWcNHDygvodSUjPYc4vEQBMH9ICB02dLSKllJubS3Z2Nl5eXgoSclkmkwkvLy+ys7PJzc2t0Lb0CWVrF/tLuF2YHrvxNWBnz0e/RJCek0+b+p4Maxtsu/pEpMa52NmyMjvYSe108TVS0Q66ChO2drG/hHHhF9m4P/Fp2czbEgnA1Oua6z8LESkXvXfIlVTWa0RhwtYunplIO2/9Gj6Aj345QUZOPm1DPLmuVaDtahMRESkFhQlbMgzrtNkAljxwD+KYpR6fbjkBwNRrdVZCRESqP4UJW0qMhOxksLMOqrE4+/Do13vIyrVwTTN/rm0VCBYL5GmMuIiIrT377LOEhYVhNpvx9vYGYMCAAQwYMOCKj42MjMRkMvHpp59WaY2zZ8+u8n0UR0NDbenoGutX74aQcJyE9Gz2Jibj5eLAa2M7YMrPhWUPQU463PoZ2Nnbtl4RkTrqu+++48UXX+SZZ55h2LBhODlZZyOePXu2jSsrbPbs2fj7+zNx4sSrul+FCVvJToWfX7V+3+x69m39gS+TBwDw0uh2BDvnwZd3wvF1YLK3XlE0rIfNyhURqcv27dsHwMMPP0xg4J992Vq3bm2rkqoVNXPYypZ3ICMOfJsQGXgtd+c8wQLLYNqGeDK8sRnmjbQGCQdXuOMrBQkRqTDDMMjIyas2N8Mwyv1cDh06xO23305QUBBOTk6EhYVx1113FUwdvm/fPm666SZ8fHxwdnamY8eOzJs3r9A2NmzYgMlk4ssvv+SZZ56hfv36eHp6ct1113H48OGC9Ro1asSzzz4LQFBQECaTiZkzZwLFN3NER0dz66234uHhgZeXF+PGjePcuXPFPo8dO3Zw44034uvri7OzM506deLrr78utM6nn36KyWRi/fr1TJ48GX9/f/z8/Bg9ejTR0dGF6ty/fz8///wzJpMJk8lEo0aNynN4y0xnJmwhLQa2vAvAj82f47Gl2aTiQ2PTWT4cGQ5zBkPiCeu1Ou5YBA262LhgEakNMnPzaf2v1bYuo8CB54fg6lj2j6E9e/bQt29f/P39ef7552nWrBlnz55l2bJl5OTkEBkZSe/evQkMDOTtt9/Gz8+Pzz//nIkTJ3L+/HmeeOKJQtt7+umn6dOnDx9//DEpKSnMmDGDG264gYMHD2Jvb8+SJUt47733mDNnDqtWrcLLy4sGDRoUW1tmZibXXXcd0dHRvPTSSzRv3pwVK1Ywbty4IuuuX7+eoUOH0qNHD95//328vLxYuHAh48aNIyMjo0hTxX333ceIESP44osviIqK4vHHH2f8+PGsW7cOgCVLljB27Fi8vLwKml8uNsdUNYUJW9j9BeSmk1evC68f9CY1N42WppO87TGf+t+8DOmx4B0G45eAf1NbVysiUq089thjmM1mfvvtNwICAgru/9vf/gbAzJkzycnJYf369YSGWi9FMHz4cJKSkpg1axb3338/Xl5eBY9r3bo1n3/+ecHP9vb23HrrrWzfvp2ePXvSqVOngvDQpUsX/P39L1vbvHnzOHjwIN999x033ngjANdffz2ZmZl89NFHhdadMmUKbdq0Yd26dZjN1o/jIUOGEBcXx9NPP81dd92Fnd2fDQhDhw7l7bffLvg5ISGBJ554gnPnzhEcHEynTp1wcXHB09OTnj17lu2gVpDChC0cXA6AufPfeMS5Gf9esIbvHP+JY44FciwQ3A7+9i14BNm4UBGpTVwc7Dnw/BBbl1HAxaHsncozMjL4+eefuffeewsFiUutW7eOa6+9tiBIXDRx4kRWrlzJ1q1bGTp0aMH9Fz/0L2rfvj0AJ0+eLPOH8vr16/Hw8CiyzTvuuKNQmDh27BiHDh3i9ddfByAvL69g2fDhw/n+++85fPgwrVq1KlWdwcG2nSlZYeJqSz4DZ3YAJmg5ku0b4hhh/ytOpgsvpEb94LYF4FzxS8KKiFzKZDKVq1mhOklMTCQ/P/+yzQwA8fHx1KtXr8j99evXL1h+KT8/v0I/X2wayMzMLHN98fHxBAUV/Ufwrx/2589bJyqcPn0606dPL3ZbcXFxVVZnZavZr6qa6ND31q9hPcEjiIGspp/DF38uH/66goSIyGX4+vpib2/P6dOnL7uOn58fZ8+eLXL/xc6KJTVTVJSfnx+//fZbkfv/2gHzYg1PPfUUo0ePLnZbLVq0qPwCq4hGc1xtF5o4aHUDAP2G3orRsDeG64XEGX/URoWJiFR/Li4u9O/fn0WLFhX5z/2ia6+9lnXr1hUa6QDw2Wef4erqWqX9CQYOHEhqairLli0rdP8XX3xR6OcWLVrQrFkz9uzZQ9euXYu9eXh4lHn/Tk5ONjlToTBxNaXHwcnN1u9bjrR+dXDGNH4JpibXWX+OOWCb2kREaog333yT3NxcevTowUcffcT69etZuHAhd9xxB6mpqTz33HM4ODgwcOBAFixYwMqVKxk/fjwrVqxg5syZhTpfVra77rqL5s2bc9ddd/Hee++xZs0apk6dyurVRUfRfPDBB6xdu5YhQ4bw5ZdfsnHjRpYuXcpLL73ELbfcUq79t2vXjj179vDVV1+xfft29u7dW9GnVCpq5riaDv8AhgXqdQCfhn/e7+AMwW1gL3B+v83KExGpCTp06MBvv/3Gc889x1NPPUVqairBwcEMGjQIR0dHWrRowZYtW3j66ad58MEHyczMpFWrVsydO7fKZ4Z0dXVl3bp1PPLIIzz55JOYTCauv/56Fi5cSO/evQutO3DgQH777TdefPFFpk6dSmJiIn5+frRu3Zpbb721XPufNWsWZ8+e5e9//zupqak0bNiQyMjISnhmJTMZFZk1pAZISUnBy8uL5ORkPD1t1BchI9F6VdDvp8KprTDoWej3eOF1jv4EC8aAf3N4aLtNyhSR2iErK4sTJ07QuHFjnJ2dbV2OVGNXeq2U9jNUZyaqSk4GbHnb2kfi/L7Cy1rdWHT9oAtTssYfh9ws69kKERGRGkBhoirEH4ev7yocIuzM0KA7tL4RAorpoetRD5y9ISsJ4g5bm0JERERqAIWJynZ4JSy+33ppcSdPyE4BsxP8Yxd4hVz+cSYTBLWxdtA8f0BhQkREagyN5qgslnxY+wJ8eZs1SDToDh4XJinpMbnkIHFR4IWmjhh1whQRkZpDZyYqQ3o8fHsvRKy3/tzjAajXEZY+AE5e0Hdq6bZzsd+ERnSIiEgNojBRUWd+h68nQHKU9XLhN74D7cbCexcmRenzMLj4lG5bgW2sX89rrgkREak5FCYqYt9iWDoZ8rLAtwmM+9x6diEvG2IPWtfpfFfptxd44YIuaecgIwFcfSu/ZhERkUqmMFEehgG/vAHrXrD+3HwojP4QnC/MqpYUZf3q4AZuxV/VrljOnjDgKfAMsY7+EBERqQH0iVVWeTnWyad2L7D+3PNBuP4FsLvkUrpJJ61fvcOsozTKYsCTlVKmiIjI1aIwUVZn98CehWCyh+GvQrf7iq6TdMr69dIps0VERGopDQ0tq9BucOPbcMfXxQcJ+DNMeIddvbpERMRmNmzYgMlkYsOGDZW2zU8//RSTyXRVrq1RUTozUR6dxpe8/NJmDhERqfU6d+7M1q1bad26ta1LsQmFiaqQGGn9qjAhIlKr5ebmYjKZ8PT0pGfPnrYux2bUzFHZLPkQc2FYaGDdTKgiUk0ZBuSkV59bBS5avWLFCjp27IiTkxONGzfm9ddfZ+bMmZgudHqPjIzEZDLx6aefFnmsyWRi5syZhe47evQod9xxB4GBgTg5OdGqVSvee++9QutcbMqYP38+06ZNIyQkBCcnJ44dO3bZZo4dO3Zw44034uvri7OzM506deLrr78uUtO2bdvo06cPzs7O1K9fn6eeeorc3NxyH5+rTWcmKltCBORmgNkFfMNtXY2IyJ9yM+D/6tu6ij89HQ2ObmV+2Nq1a7npppvo1asXCxcuJD8/n1dffZXz58+Xq4wDBw7Qu3dvwsLCeOONNwgODmb16tU8/PDDxMXF8dxzzxVa/6mnnqJXr168//772NnZERgYyLlz54psd/369QwdOpQePXrw/vvv4+XlxcKFCxk3bhwZGRlMnDixYP/XXnstjRo14tNPP8XV1ZXZs2fzxRdflOv52ILCRGU7t9f6Nah14eGiIiJSKZ555hmCgoL48ccfcXZ2BmDIkCE0atSoXNt77LHH8PDwYNOmTXh6egIwePBgsrOzefnll3n44Yfx8flzJuMmTZqwaNGiK253ypQptGnThnXr1mE2mwvqjIuL4+mnn+auu+7Czs6O559/HsMwWLduHUFBQQCMGDGCtm3bluv52ILCRGW7GCaC29m2DhGRv3JwtZ4NqC4cXMv8kPT0dLZv386UKVMKggSAh4cHN9xwA/PmzSvT9rKysli7di2TJ0/G1dWVvLy8gmXDhw/n3XffZdu2bQwbNqzg/jFjxlxxu8eOHePQoUO8/vrrAEW2+/3333P48GFatWrF+vXrufbaawuCBIC9vT3jxo1j1qxZZXo+tqIwUdnO77N+Dao5iVJE6giTqVzNCtVJYmIiFouF4ODgIsuKu+9K4uPjycvL45133uGdd94pdp24uLhCP9erV++K273Y5DJ9+nSmT59e4nbj4+Mr7fnYisJEZTt3IUwEt7dtHSIitZCPjw8mk6nYPgqX3nfxrEV2dnahdeLj44tsz97enjvvvJMHH3yw2H02bty40M+mUsxs7O/vD1j7V4wePbrYdVq0aAGAn5/fFZ9PdacwUZnS4yH1winEII3kEBGpbG5ubnTv3p3Fixfz2muvFYSG1NRUli9fXrBeUFAQzs7O/PHHH4Ue/9133xX62dXVlYEDB7Jr1y7at2+Po6NjpdTZokULmjVrxp49e/i///u/EtcdOHAgy5Yt4/z58wVNHfn5+Xz11VeVUsvVoDBRmc5f6C/hGw5OHratRUSklnrhhRcYOnQogwcPZtq0aeTn5/PKK6/g5uZGQkICYD17MH78eD755BOaNGlChw4d+O2334odIfHf//6Xvn37cs011zB58mQaNWpEamoqx44dY/ny5axbt65cdX7wwQcMGzaMIUOGMHHiREJCQkhISODgwYPs3LmzoBPns88+y7Jlyxg0aBD/+te/cHV15b333iM9Pb38B+kqU5ioTAUjOdRfQkSkqgwePJilS5fy7LPPMm7cOIKDg5kyZQqZmZmFOiy+8cYbALz66qukpaUxaNAgvv/++yKjPlq3bs3OnTt54YUXePbZZ4mJicHb25tmzZoxfPjwctc5cOBAfvvtN1588UWmTp1KYmIifn5+tG7dmltvvbVgvbZt2/LTTz8xbdo0JkyYgI+PD3feeSdjxoxh0qRJ5d7/1WQyjArMGlIDpKSk4OXlRXJycsGQnyqTkw7nD4C9Gep3qtp9iYhcRlZWFidOnKBx48aFRjzUdjNnzmTWrFnU8o+1SnWl10ppP0N1ZqIyObpZLwQmIiJSh2g6bREREakQhQkREakVZs6cqSYOG1GYEBERkQpRmBAREZEKUZgQEamldMpfrqSyXiMKEyIitYyDgwMmk6lGTXoktpGeno7JZMLBwaFC29HQUBGRWsbe3h4vLy9iY2PJzs7G09MTs9lcqmtKSO1nGAZ5eXmkpKSQkpKCt7c39vb2FdqmwoSISC0UHByMi4sLMTExpKSk2LocqYbs7e2pV68eXl5eFd6WwoSISC1kMpnw9vbGy8uL/Px88vLybF2SVCNmsxl7e/tKO1ulMCEiUouZTCbMZjNms97upeqoA6aIiIhUiMKEiIiIVIjChIiIiFSIwoSIiIhUSK3vkXNxdi8NjRIRESmbi5+dV5ops9aHidTUVABCQ0NtXImIiEjNlJqaWuJ8FCajlk/ebrFYiI6OxsPDQ7O/VaKUlBRCQ0OJiorC09PT1uXUKjq2VUfHturo2FYdWx5bwzBITU2lfv362NldvmdErT8zYWdnR4MGDWxdRq3l6empN44qomNbdXRsq46ObdWx1bEtzQyZ6oApIiIiFaIwISIiIhWiMCHl4uTkxHPPPYeTk5OtS6l1dGyrjo5t1dGxrTo14djW+g6YIiIiUrV0ZkJEREQqRGFCREREKkRhQkRERCpEYUJEREQqRGGijtq9ezcjRowgLCwMFxcXfH196dWrF59//nmRdXNzc3nzzTdp164dLi4ueHt707t3b7Zs2XLF/QwYMACTyVTkNnTo0Kp4WtVCaY9tccfl4q1ly5al2tdPP/1Er169cHV1xd/fn4kTJxITE1MVT6tauFrHVq/byx9bwzD46KOP6NKlC56envj5+dG/f39WrFhR6n3pdVs1x9aWr9taPwOmFC8pKYnQ0FBuv/12QkJCSE9PZ8GCBdx5551ERkby7LPPApCfn8+oUaPYtGkTTzzxBL179yY9PZ3ff/+d9PT0Uu0rPDycBQsWFLrP29u7sp9StVHaY7t169Yij/3111+ZOnUqo0aNuuJ+fv75Z4YNG8aIESP47rvviImJYcaMGVx77bXs2LGjWg8jK6+rdWxBr9vLHdvnnnuOF154gQceeICXX36ZrKws3nnnHUaOHMm3337L6NGjS9yPXrdVd2zBhq9bQ+QSPXr0MEJDQwt+/s9//mPY2dkZW7duLdf2+vfvb7Rp06ayyqvR/npsizNx4kTDZDIZR48eveL2unXrZrRu3drIzc0tuG/z5s0GYMyePbvC9dYklX1s9br901+PbUhIiNG3b99C62RmZhpeXl7GjTfeeMXt6XX7p8o+trZ83aqZQwrx9/fHbP7zhNV///tf+vXrR8+ePW1YVe3w12P7V6mpqSxatIj+/fvTtGnTErd15swZtm/fzp133llom71796Z58+YsWbKk0uquCSrz2Ephfz22Dg4ORa7V4OzsXHAriV63hVXmsbU1hYk6zmKxkJeXR2xsLLNnz2b16tXMmDEDgKioKCIjI2nXrh1PP/00QUFBmM1m2rRpw7x580q9j+PHj+Pr64vZbKZJkyY888wzZGZmVtVTqjZKOrbFWbhwIenp6dx3331X3Pa+ffsAaN++fZFl7du3L1heW1Xlsb1Ir9vij+0jjzzCqlWrmDNnDomJiZw9e5bHHnuM5ORkHn744RK3rddt1R3bi2z2urXJ+RCpNu6//34DMADD0dGx0GnGrVu3GoDh6elptG7d2vj666+N1atXG2PHjjUA48MPP7zi9p955hlj9uzZxrp164wVK1YYDz30kGE2m41+/foZ+fn5VfnUbK6kY1ucHj16GN7e3kZmZuYVt71gwQIDKLb5adKkSYajo2O5664JqvLYGoZet1c6tu+//77h5ORUsJ6vr6/x448/XnHbet1W3bE1DNu+bhUm6riTJ08a27dvN1asWGE88MADhp2dnfHaa68ZhvFnO6ajo6MRGRlZ8BiLxWJ07tzZaNCgQbn2+frrrxuAsXjx4kp5DtVVScf2r/bt22cAxoMPPliqbV98U962bVuRZZMmTTKcnJwqVHt1V5XH9nL0urX65JNPDCcnJ2PatGnGTz/9ZPzwww/GbbfdZri6uhqrVq0qcdt63Vbdsb2cq/W6VZiQQh544AHDbDYbMTExxqFDhwzAaN++fZH1nnrqKQMwzp8/X+Z9nDt3zgCMJ554ojJKrjEuPbZ/9eijjxqAsWvXrlJta9WqVQZgrFixosiysWPHGvXq1atouTVKZR7by9HrNsZISEgwXFxcig1m/fv3Nxo1alTitvS6Lawyj+3lXK3XrfpMSCHdu3cnLy+PiIgImjRpgqura7HrGReuD2dnV/6XUEUeWxNdemwvlZOTw/z58+nSpQsdO3Ys1bbatm0LwN69e4ss27t3b8HyuqIyj+2V1OXX7eHDh8nMzKRbt25F1uvatSuRkZGkpaVddlt63RZWmcf2Sqr6dVu3/irkitavX4+dnR3h4eGYzWZuuukmDh48SGRkZME6hmGwatUqmjRpgr+/f5n3cbHzZl0bIXLpsb3UsmXLiIuL49577y31tkJCQujevTuff/45+fn5Bfdv27aNw4cPl2o8em1Smcf2cvS6Dad+/fqA9XV2KcMw2LZtGz4+Pri5uV12W3rdFlaZx/ZyrtrrtkrPe0i19fe//92YNm2a8dVXXxkbNmwwvvnmG2PcuHEGYDz++OMF6x07dszw9vY2WrRoYXz55ZfGihUrjFGjRhkmk8lYtGhRoW3a29sbgwYNKvh548aNxpAhQ4z333/fWLNmjbFs2TJj8uTJBevV1o5spT22Fw0dOtRwcXExkpKSLrvNvx5bwzCM9evXG2az2Rg1apTx448/GgsWLDBCQ0ONtm3bGllZWZX+vKqDq3Fs9bot+diOHj3asLOzMx555BFj9erVxrJly4wxY8YYgPHCCy8U2qZet1ZX49ja+nWrMFFHffLJJ8Y111xj+Pv7G2az2fD29jb69+9vzJ8/v8i6e/fuNUaMGGF4eHgYzs7ORs+ePY3ly5cXWQ8w+vfvX/Dz0aNHjeHDhxshISGGk5OT4ezsbLRr18548cUXa+2bhmGU7dieOnXKsLOzM+66664St/nXY3vRmjVrjJ49exrOzs6Gr6+vcdddd5WrH0tNcTWOrV63JR/bzMxM47XXXjPat29veHh4GL6+vkbPnj2Nzz//3LBYLIXW1evW6mocW1u/bk0XihIREREpF/WZEBERkQpRmBAREZEKUZgQERGRClGYEBERkQpRmBAREZEKUZgQERGRClGYEBERkQpRmBAREZEKUZgQERGRClGYEBERkQpRmBAREZEK+X8smx1m1icWHQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 600x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot\n",
    "flag1, flag2 = 0, 0\n",
    "fig, ax = plt.subplots(figsize=[6, 6])\n",
    "for line, line_id in zip(line_strs, line_ids):\n",
    "    kwargs = dict(color = 'C0' if line_id == 1 else 'C1')\n",
    "    # label the first fault segment for legend\n",
    "    if not flag1 and line_id == 1:\n",
    "        kwargs['label'] = 'confident'\n",
    "        flag1 = 1\n",
    "    elif not flag2 and line_id == 2:\n",
    "        kwargs['label'] = 'queried'\n",
    "        flag2 = 1\n",
    "    ax.plot(line.coords.xy[0], line.coords.xy[1], **kwargs)\n",
    "# axis format\n",
    "ax.set_title('2023 Turkey EQ Simple Fault')\n",
    "ax.set_aspect('equal')\n",
    "ax.legend()\n",
    "# save\n",
    "out_fig = os.path.abspath(f'../fault_{version}.png')\n",
    "print('save figure to file', out_fig)\n",
    "plt.savefig(out_fig, bbox_inches='tight', transparent=True, dpi=300)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2e7b428a-af0d-493c-abd3-47194160683b",
   "metadata": {},
   "source": [
    "### Save to GMT lonlat file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "8f27436f-3796-4d2b-88c8-32f0cb15809c",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "write file: /Users/yunjunz/data/archives/2023TurkeyEQ/USGS/simple_fault_2023-03-15/../simple_fault_confident.lonlat\n",
      "write file: /Users/yunjunz/data/archives/2023TurkeyEQ/USGS/simple_fault_2023-03-15/../simple_fault_queried.lonlat\n"
     ]
    }
   ],
   "source": [
    "for suffix, ref_id in zip(['_confident', '_queried'], [1, 2]):\n",
    "    out_file = os.path.join(work_dir, f'../simple_fault{suffix}.lonlat')\n",
    "    num_line = 0\n",
    "    with open(out_file, 'w') as f:\n",
    "        for line_str, line_id in zip(line_strs, line_ids):\n",
    "            if line_id == ref_id:\n",
    "                # new line - start\n",
    "                num_line += 1\n",
    "                f.write(f'> segment {num_line}\\n')\n",
    "                xs, ys = line_str.coords.xy\n",
    "                for x, y in zip(xs, ys):\n",
    "                    f.write(f'{x:.9f}\\t{y:.9f}\\n')\n",
    "    print(f'write file: {out_file}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ebe10532-a629-4b82-8e7f-2fe72239c785",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
